2.1 Programme

Wir brauchen zwei Programme:

1. R
   
2. RStudio (Desktop-Version)

Bitte laden Sie diese herunter und installieren Sie sie. Wenn R installiert ist, dann findet RStudio R auch direkt.

Zur Installation gehen Sie so vor:

1. Laden Sie R herunter und öffnen Sie die heruntergeladene Installationsdatei; folgen Sie den Hinweisen, die Sie durch die Installation leiten
   • Windows hier (Windows 7 oder neuer)
   • Mac hier (OSX 10.11 oder neuer) ⚠️ wählen Sie die neueste R-Version (höchste Versionsnummer)
   • Linux hier
2. Laden Sie RStudio herunter (Desktop-Version)
   • Sie finden eine Version für alle gängigen Betriebssysteme.

Beide Programme sind kostenlos.

Wenn alles läuft, sieht es etwa so aus:

2.2 Hilfe

Bei Sauer (2019) findet sich ein Abschnitt, der bei gängigen R-Problemen weiterhilft (Abschnitt 3.8; s. auch 3.1).

Hier finden sich einige Einführungen in R in unterschiedlichem Niveau; Antworten auf häufige Fragen finden sich hier. Außerdem hilft erfahrungsgemäß: Googeln Sie nach der Fehlermeldung.

2.3 Erweiterungen (Pakete, engl. packages)

Um einen Befehl zu verwenden, der nicht im Standard-R, sondern in einer Erweiterung von R (“Paket”) wohnt, müssen sie dieses Paket erst starten (laden). Dazu können Sie den Befehl library verwenden. Wir benötigen diese Pakete; bitte laden (d.h. Code ausführen); starten Sie jetzt diese Erweiterungen (per Klick oder mit folgenden Befehlen).

library(mosaic)  # Zugpferd
library(openxlsx)  # Excel-Dateien schreiben
library(GGally)  # Korrelationsdiagramme
library(lsr)  # Effektstärkemaße 
library(sjmisc)  # Datenjudo
library(tidyverse)  # Datenjudo

Falls Sie R mit einem beleidigten Kommentar konfrontiert wie “could not find function XXX”, dann liegt es vermutlich daran, dass Sie vergesse haben, das Paket, in dem die Funktion XXX “wohnt”, zu laden.

Oder Sie klicken den Namen des Pakets hier an:

Wir gehen im Folgenden davon aus, dass Sie diese beiden Pakete geladen haben.

Nach jedem Start von R bzw. RStudio müssen Sie die Erweiterung erneut laden (wenn Sie sie benutzen wollen).

⚠️ Um ein Paket zu laden, muss es installiert sein. Klicken Sie zum Installieren auf den Button “Install” unter dem Reiter “Packages” in RStudio:

Sie müssen ein Paket nur einmal installieren, um es verwenden zu können. Sie installieren ja auch nicht Ihren Browser jedes Mal neu, wenn Sie den Computer starten.

2.4 Daten

Wir verwenden in diesem Kurs diese Datensätze:

• TeachingRatings; Sie können ihn hier herunterladen (XLSX-Version: https://data-se.netlify.com/download/TR.xlsx).
• mtcars; mtcars ist schon im Standard-R fest eingebaut; Sie müssen also nichts weiter tun.
• tips; den Datensatz tips können Sie hier herunterladen.

⚠️ Bitte stellen Sie sicher, dass Sie auf diese Daten zugreifen können während des Kurs. Laden Sie sie vorab herunter.

4.1 Kringel-Notation (Tilde; Formel-Schreibweise)

Das Paket mosaic wird unser Zugpferd für alle folgenden Analysen ein. Es hat den Charme, über eine einfache, konsistente Syntax zu verfügen. Mit wenig kann man da viel erreichen. Genau das, wovon man als Student träumt… (so denken Dozenten, jedenfalls).

Die folgende Syntax

Zielbefehl(y ~ x , data=...)

wird verwendet für

• graphische Zusammenfassungen,
• numerische Zusammenfassungen und
• inferentstatistische Auswertungen

Für Grafiken gilt:

• y: y-Achse Variable
• x: x-Achse Variable

Generell gilt:

y ~ x

DAs kann in der Regel gelesen werden y hängt ab von x. Die “Kringel-Schreibweise” nennt man in R auch eine “Formel” (formula) oder entsprechend das “formula interface”. Der Kringel (die Tilde) ~ erzeugt sich beim Mac mit ALT+n und bei Windows steht es auf einer Taste ziemlich rechts auf der Tastatur. Die Verwendung der Tilde wird auch als “Formel-Schreibweise” (engl. “Formula Interface”) bezeichnet.

5.1 tidy data - Tabellen in Normalform

Damit Sie in R vernünftig mit Ihren Daten arbeiten können, sollten die Daten “tidy” sein, d.h. in Normalform. Was ist Normalform? Betrachten Sie folgende Abbildung - so sieht eine Tabelle in Normalform aus.

Übrigens heißen Tabellen (mit Spaltennamen) in R Dataframes.

Die goldene Regel der Normalform einer Tabelle lautet also:

In jeder Zeile steht eine Beobachtung (z.B. Person). In jeder Spalte eine Variable (z.B. Geschlecht). In der ersten Zeile stehen die Spaltennamen, danach folgen die Werte. Sonst steht nichts in der Tabelle.

⚠️ Falls Ihre Daten nicht in Normalform sind, sollten Sie diese zunächst in Normalform bringen.

💡 Der einfachste Weg (von der Lernkurve her betrachtet, nicht vom Zeitaufwand), Daten in Normalform zu bringen, ist sie in Excel passend umzubauen.

5.2 Beispiel für Daten in Nicht-Normalform

Sie denken, dass Ihre Daten immer/auf jeden Fall in Normalform sind? Dann schauen Sie sich mal dieses Bild an:

Wir werden in diesem Kurs nicht bearbeiten, wie man Daten von “breit” auf “lang” (=tidy) umformatiert. Aber lesen Sie bei Interesse doch z.B. hier nach.

5.3 Daten anschauen

Es empfiehlt sich, zu Beginn einen Blick auf die Daten zu werfen, um zu prüfen, ob alles augenscheinlich seine Richtigkeit hat. Tun Sie das immer, viel Ärger lässt sich so ersparen.

glimpse(TeachingRatings)
#> Rows: 463
#> Columns: 13
#> $ X           <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
#> $ minority    <fct> yes, no, no, no, no, no, no, no, no, no, yes, no, no, no,…
#> $ age         <int> 36, 59, 51, 40, 31, 62, 33, 51, 33, 47, 35, 37, 42, 49, 3…
#> $ gender      <fct> female, male, male, female, female, male, female, female,…
#> $ credits     <fct> more, more, more, more, more, more, more, more, more, mor…
#> $ beauty      <dbl> 0.2899, -0.7377, -0.5720, -0.6780, 1.5098, 0.5886, -0.126…
#> $ eval        <dbl> 4.3, 4.5, 3.7, 4.3, 4.4, 4.2, 4.0, 3.4, 4.5, 3.9, 3.1, 4.…
#> $ division    <fct> upper, upper, upper, upper, upper, upper, upper, upper, u…
#> $ native      <fct> yes, yes, yes, yes, yes, yes, yes, yes, yes, yes, no, yes…
#> $ tenure      <fct> yes, yes, yes, yes, yes, yes, yes, yes, yes, no, yes, no,…
#> $ students    <int> 24, 17, 55, 40, 42, 182, 33, 25, 48, 16, 18, 30, 28, 30, …
#> $ allstudents <int> 43, 20, 55, 46, 48, 282, 41, 41, 60, 19, 25, 34, 40, 36, …
#> $ prof        <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…

5.4 Selber Daten erzeugen

Normalerweise werden Sie Daten in R einlesen, aber manchmal möchte man selber Daten erzeugen. Dazu sollten Sie wissen: Ein Dataframe besteht aus Spalten, wie Sie wissen. Diese Spalten nennt man auch Vektoren. Ein Vektor ist eine Sammlung von einzelnen Datenstücken, die hintereinander aufgereiht sind, wie Wäsche an der Wäscheleine. Ach ja, Vektoren müssen “sortenrein” sein, also nur Zahlen oder nur Text etc. Um einen Vektor zu erzeugen, benutzt man den Befehle c (wie combine oder concatenate, aneinanderfügen):

meine_freunde <- c("Bibi", "Uschi", "Ulf", "Donald der Große")

Deine_PINs <- c(1234, 7777, 1234567, 0000)

💻 AUFGABE:

1. Erzeugen Sie einen Vektor mit einem kreativen Namen.
2. Erzeugen Sie einen Vektor, in dem Sie Zahlen und Text mischen. Was passiert?
3. Erstellen Sie einen Vektor aus meine_freunde und Deine_PINs. Was passiert?
4. Adeln Sie mal Ihren Vektor zu einer waschechten Tabelle (Dataframe) mit dem Befehl data.frame(vektor); nicht vergessen, der resultierenden Tabelle einen Namen zu geben bzw. eine neue Tabelle zu benennen, die das Ergebnis des Befehle data.frame speichert.

5.5 Daten als CSV oder Excel exportieren

Wie kriege ich die Daten aus R wieder raus? Was ist sozusagen mit REXIT-Strategie? Keine Sorge, die Straße fährt in beide Richtungen. Sagen wir, Sie möchten den Dataframe TeachingRatings exportieren, um außerhalb von R damit Kunststücke zu vollbringen:

write.csv(TeachingRatings, "TeachingRatings.csv")

write.xlsx(TeachingRatings, "TeachingRatings.xlsx")

Dieser Befehl (write.xlsx) schreibt eine XLSX-Datei in das aktuelle R-Verzeichnis (das Arbeitsverzeichnis).

5.6 Das R-Arbeitsverzeichnis

💡 R schreibt Dateien immer in das R-Arbeitsverzeichnis. Das R-Arbeitsverzeichnis ist das Arbeitsverzeichnis, das R als Standard annimmt. Wenn Sie zu R etwas sagen wie “lade daten.csv!”, dann wird R in dem Ordner nach schauen, der als Arbeitsverzeichnis ausgewählt ist.

Mit getwd() findet man heraus, was das aktuelle Arbeitsverzeichnis ist:

getwd()
#> [1] "/Users/sebastiansaueruser/Documents/Lehre/Statistik/r-crashkurse"

Mit set.wd() kann man ein Arbeitsverzeichnis wählen. Alternativ geht das auch über das Menü Session > Set Working Directory > Choose Directory ….

💻 AUFGABE:

1. Finden Sie Ihr aktuelles Arbeitsverzeichnis heraus.
2. Setzen Sie Ihr Arbeitsverzeichnis auf den Ordner, in dem Ihre Lieblings-Skriptdatei liegt.

5.7 Textkodierung in UTF-8

Falls Sie RStudio oder ein beliebiger Texteditor irgendwann fragt, wie die Textdatei kodiert sein soll, wählen Sie immer “UTF-8”. UTF-8 ist eine Kodierungstabelle, so dass der Computer weiß, welcher Buchstabe welcher Bitfolge zugeordnet ist; dabei ist UTF-8 so großzügig geplant, dass alle möglichen Sonderzeichen (Deutsch, Chinesisch, Hebräisch…) dazu passen. UTF-8 ist die Standard-Kodierung für Textdateien im Internet. Falls R-Studio Sie nach “Save with Encoding” fragt, wählen Sie immer UTF-8. Speichern Sie Ihre Skriptdateien am besten immer in diesem Format. Etwas mehr zur Textkodierung finden Sie bei Sauer (2019) oder hier.

💡 In RStudio kann man unter File..Save with Encoding... die Textcodierung einstellen.

💻 AUFGABE:

1. Prüfen Sie, in welchem Format Ihr Dokument kodiert ist.
2. Setzen Sie das Format ggf. auf UTF-8.

6.1 Auf Fehler prüfen beim Einlesen

⚠️ Ein häufiger Fehler ist, dass die Daten nicht richtig eingelesen werden. Zum Beispiel werden die Spaltentrennzeichen nicht richtig erkannt. Das kann dann so aussehen:

Unter “delimiter” in der Maske können Sie das Trennzeichen anpassen.

⚠️ “Deutsche” CSV-Dateien verwenden als Dezimaltrennzeichen ein Komma; englisch-formatierte CSV-Dateien hingegen einen Punk. R geht per Default von englisch-formatierten CSV-Dateien aus. Importieren Sie eine deutsch-formatierte CSV-Datei, müssen Sie das Dezimaltrennzeichen von Hand ändern; es wird nicht automatisch erkannt.

💡 Unter “locale” können Sie das Dezimaltrennzeichen ggf. anpassen.

6.2 Spaltennamen korrigieren

Spaltennamen müssen auch “tidy” sein. Das heißt in diesem Fall:

• keine Leerzeichen
• keine Sonderzeichen (#,ß,ä,…)
• nicht zu lang, aber trotzdem informativ

Spaltennamen sollten nur Buchstaben (ohne Umlaute) und Ziffern enthalten und mit Buchstaben beginnen. Für Textdaten in den Spalten sind diese Regeln auch sinnvoll.

💡 Am einfachsten ändern Sie die Spaltennamen in Excel.

In R können Sie Spaltennamen z.B. so ändern:

rename(TeachingRatings, festangestellt = tenure) -> TR2

In Pseudo-R könnte man schreiben:

benenne_spalte_um(meine_tabelle, 
                  neuer_name = altername) -> 
  meine_neue_tabelle

Der R-Zufweisungspfeil <- bzw. -> funktioniert in beide Richtungen; er darf nach links oder rechts zeigen. In jedem Fall wird das Objekt, auf das er zeigt, “befüllt” mit den Inhalten die auf der anderen Seite stehen.

💻 AUFGABE:

• Benennen Sie in TeachingRatings die Spalte native in Muttersprachler um; speichern Sie aber das Ergebnis in einem neuen Dataframe.
• Suchen Sie sich noch zwei weitere Spalten, und benennen Sie die Spaltennamen nach eigenen Vorstellungen um!

6.3 Umkodieren

Gerade bei der Analyse von Fragebogendaten ist es immer wieder nötig, Daten umzukodieren. Klassisches Beispiel: Ein Item ist negativ kodiert. Zum Beispiel das Item “Ich bin ein Couch-Potator” in einem Fragebogen für Extraversion.

Nehmen wir an, das Item “i04” hat die Werte 1 (“stimme überhaupt nicht zu”) bis 4 (“stimme voll und ganz zu”). Kreuzt jemand das Couch-Potato-Item mit 4 an, so sollte er nicht die maximale Extraversion-Punktzahl (4), sondern die minimale Extraversion-Punktzahl (1) erhalten. Also

` 1 –> 4

2 –> 3

3 –> 2

4 –> 1 `

Am einfachsten ist dies zu bewerkstelligen mit folgendem R-Befehl aus dem Paket sjmisc.

Zuerst erzeugen wir uns ganzzahlige Evaluationswerte:

TeachingRatings <- mutate(TeachingRatings, eval_rounded = round(eval))
tally(~ eval_rounded, data = TeachingRatings)
#> eval_rounded
#>   2   3   4   5 
#>   4  73 307  79

Das war die Vorbereitung; jetzt das eigentliche Umkodieren:

TeachingRatings <- rec(TeachingRatings, eval_rounded, rec = "2=5; 3=4; 4=3; 5=2")

rec() kodiert die angegebene Spalte, hier eval_rounded um, genau so wie es mit rec beschrieben ist. Dabei wird eine neue Spalte (Variable) angelegt, deren Name der umzukodierenden Variable plus dem Suffix _r entspricht, hier also eval_rounded_r (s. Tabelle ??)

6.4 Fehlende Werte

Der einfachste Umgang mit fehlenden Werten ist: nichts machen. Denken Sie nur daran, dass viele R-Befehle von Natur aus nervös sind - beim Anblick von fehlenden Werten werden sie panisch und machen nix mehr. Strecken alle Viere von sich. Der einfachste Weg ist, Zeilen mit fehlenden Werten (in den betreffenden Spalten) zu entfernen. Das geht zum Beispiel mit drop_na() aus tidyverse²:

mein_bereinigter_dataframe <- drop_na(mein_dataframe,
                                      spalte1, spalte2)

Dieser Befehl löscht alle Zeilen mit fehlenden Werten aus der Tabelle³ mein_dataframe, sofern sie sich in den Spalten spalte1 und spalte2 befinden (Sie können beliebig viele Spalten anführen; einfach mit Komma trennen. Geben Sie keine Spalten an, so werden alle Spalten nach fehlenden Werten durchsucht). Das Ergebnis ist eine Tabelle ohne fehlende Werte in den Spalten spalte1 und spalte2; wir speichern es als neue Tabelle mit Namen mein_bereinigter_dataframe.

Alternativ kann man einzelnen Befehlen, wie mean() sagen, dass sie fehlende Werte ignorieren sollen. Haben sie fehlende Werte in ihrer Tabelle, so verwenden Sie den Parameter na.rm = TRUE. na steht für “not available”, also fehlende Werte. rm steht für “remove”. Also mean(~ i04_r, data = meine_tabelle, na.rm = TRUE). Leider ist die genaue Schreibweise nicht immer die gleiche; hier die wichtigsten Befehle mit dem Parameter, der sei auch bei fehlenden Werten wieder zurück zur Arbeit bringt:

mean(~ eval, data = TeachingRatings, 
     na.rm = TRUE) 
#> [1] 4
# für sd, median etc. genauso

cor(eval ~ beauty, data = TeachingRatings, 
    use = "complete.obs")
#> [1] 0.189

💡 Der R-Befehl inspect aus mosaic zeigt Ihnen an, ob es fehlende Werte gibt: inspect(meine_daten).

💻 AUFGABE:

• Prüfen Sie, ob es im Datensatz TeachingRatings fehlende Werte gibt.
• Prüfen Sie, ob es im Datensatz mtcars fehlende Werte gibt.

6.5 Komische Werte

Hat ein Spaßvogel beim Alter 999 oder -1 angegeben, kann das Ihre Daten ganz schön verhageln. Prüfen Sie die Daten auf komische Werte. Der einfachste Weg ist, sich die Daten in Excel anzuschauen. Cleverer ist noch, sich Zusammenfassungen auszugeben, wie der kleinste oder der größte Wert, oder der Mittelwert etc., und dann zu schauen, ob einem etwas spanisch vorkommt. Diagramme sind ebenfalls hilfreich. Dann ändern Sie die Werte in Excel und laden die Daten erneut ins R.

6.6 Logische Variablen bilden

Sagen wir, uns interessiert welches Auto mehr als 200 PS hat; wir wollen Autos mit mehr als 200 PS vergleichen (“Spass”) mit schwach motorisierten Autos (“Kruecke”). Wie können wir das (einfach) in R erreichen? Logische Variablen sind ein einfacher Weg.

TeachingRatings <- TeachingRatings %>%
  mutate(Traumdozent = beauty > 1)

Dieser Befehl hat eine Spalte (Variable) in der Tabelle TeachingRatings erzeugt, in der TRUE steht, wenn das Auto der jeweiligen Spalte die Bedingung (beauty > 1) erfüllt⁴. Nicht glauben, nachschauen.

prop(~Traumdozent, data = TeachingRatings)
#> prop_TRUE 
#>     0.145

Ok, etwa 15% der Dozenten sind so hübsch. Erzeugen wir einen Teil-Datensatz nur mit diesen Dozentenmodells:

Dozimodels <- filter(TeachingRatings, Traumdozent == TRUE)
glimpse(Dozimodels)
inspect(Dozimodels)

💻 AUFGABE:

• Erstellen Sie eine Variable Asbach, definiert als TRUE, wenn age < 70.
• Erstellen Sie eine Variable keiner_mag_mich, definiert als TRUE, wenn eval <= 2.
• Denken Sie sich noch selber mindestens ein Beispiel aus.

6.7 Daten univariat zusammenfassen

Deskriptive Statistik ist letztlich nichts anderes, als Daten geschickt zusammenzufassen. Praktisch wird meistens eine Spalte einer Tabelle zu einer Zahl zusammengefasst.

Schauen wir uns das mal mit echten Daten an. Der Datensatz TeachingRatings ist schon in R eingebaut, so dass wir in nicht extra laden müssen. Ganz praktisch. Dazu fragen wir den Inspektor inspect, der würde uns auch noch verraten wie die nominalen Variablen sich so verteilen – wenn wir hier welche hätten.

inspect(TeachingRatings)
#> 
#> categorical variables:  
#>          name   class levels   n missing
#> 1    minority  factor      2 463       0
#> 2      gender  factor      2 463       0
#> 3     credits  factor      2 463       0
#> 4    division  factor      2 463       0
#> 5      native  factor      2 463       0
#> 6      tenure  factor      2 463       0
#> 7 Traumdozent logical      2 463       0
#>                                    distribution
#> 1 no (86.2%), yes (13.8%)                      
#> 2 male (57.9%), female (42.1%)                 
#> 3 more (94.2%), single (5.8%)                  
#> 4 upper (66.1%), lower (33.9%)                 
#> 5 yes (94%), no (6%)                           
#> 6 yes (78%), no (22%)                          
#> 7 FALSE (85.5%), TRUE (14.5%)                  
#> 
#> quantitative variables:  
#>             name   class   min      Q1  median      Q3    max     mean      sd
#> 1              X integer  1.00 116.500 232.000 347.500 463.00 2.32e+02 133.801
#> 2            age integer 29.00  42.000  48.000  57.000  73.00 4.84e+01   9.803
#> 3         beauty numeric -1.45  -0.656  -0.068   0.546   1.97 6.26e-08   0.789
#> 4           eval numeric  2.10   3.600   4.000   4.400   5.00 4.00e+00   0.555
#> 5       students integer  5.00  15.000  23.000  40.000 380.00 3.66e+01  45.018
#> 6    allstudents integer  8.00  19.000  29.000  60.000 581.00 5.52e+01  75.073
#> 7           prof integer  1.00  20.000  44.000  70.500  94.00 4.54e+01  27.509
#> 8   eval_rounded numeric  2.00   4.000   4.000   4.000   5.00 4.00e+00   0.603
#> 9 eval_rounded_r numeric  2.00   3.000   3.000   3.000   5.00 3.00e+00   0.603
#>     n missing
#> 1 463       0
#> 2 463       0
#> 3 463       0
#> 4 463       0
#> 5 463       0
#> 6 463       0
#> 7 463       0
#> 8 463       0
#> 9 463       0

💡 mit help(Befehl) bekommt man Hilfe zu einem Befehl oder einem sonstigen Objekt (z.B. Datensatz).

favstats(~ eval, data = TeachingRatings)
#>  min  Q1 median  Q3 max mean    sd   n missing
#>  2.1 3.6      4 4.4   5    4 0.555 463       0

favstats(eval ~ gender, data = TeachingRatings)
#>   gender min  Q1 median  Q3 max mean    sd   n missing
#> 1 female 2.3 3.6   3.90 4.3 4.9 3.90 0.539 195       0
#> 2   male 2.1 3.7   4.15 4.5 5.0 4.07 0.557 268       0

mean(eval~gender, data = TeachingRatings)
#> female   male 
#>   3.90   4.07
median(eval~gender, data = TeachingRatings)
#> female   male 
#>   3.90   4.15
sd(eval~gender, data = TeachingRatings)
#> female   male 
#>  0.539  0.557
var(eval~gender, data = TeachingRatings)
#> female   male 
#>   0.29   0.31
IQR(eval~gender, data = TeachingRatings)
#> female   male 
#>    0.7    0.8
diffmean(eval~gender, data = TeachingRatings)
#> diffmean 
#>    0.168
min(eval~gender, data = TeachingRatings)
#> female   male 
#>    2.3    2.1
max(eval~gender, data = TeachingRatings)
#> female   male 
#>    4.9    5.0

diff(median(eval ~ gender, data = TeachingRatings))
#> male 
#> 0.25

mean(eval~gender, data = TeachingRatings, na.rm = TRUE)
#> female   male 
#>   3.90   4.07

tally(~gender, data = TeachingRatings)
#> gender
#> female   male 
#>    195    268

inspect(TeachingRatings)
#> 
#> categorical variables:  
#>          name   class levels   n missing
#> 1    minority  factor      2 463       0
#> 2      gender  factor      2 463       0
#> 3     credits  factor      2 463       0
#> 4    division  factor      2 463       0
#> 5      native  factor      2 463       0
#> 6      tenure  factor      2 463       0
#> 7 Traumdozent logical      2 463       0
#>                                    distribution
#> 1 no (86.2%), yes (13.8%)                      
#> 2 male (57.9%), female (42.1%)                 
#> 3 more (94.2%), single (5.8%)                  
#> 4 upper (66.1%), lower (33.9%)                 
#> 5 yes (94%), no (6%)                           
#> 6 yes (78%), no (22%)                          
#> 7 FALSE (85.5%), TRUE (14.5%)                  
#> 
#> quantitative variables:  
#>             name   class   min      Q1  median      Q3    max     mean      sd
#> 1              X integer  1.00 116.500 232.000 347.500 463.00 2.32e+02 133.801
#> 2            age integer 29.00  42.000  48.000  57.000  73.00 4.84e+01   9.803
#> 3         beauty numeric -1.45  -0.656  -0.068   0.546   1.97 6.26e-08   0.789
#> 4           eval numeric  2.10   3.600   4.000   4.400   5.00 4.00e+00   0.555
#> 5       students integer  5.00  15.000  23.000  40.000 380.00 3.66e+01  45.018
#> 6    allstudents integer  8.00  19.000  29.000  60.000 581.00 5.52e+01  75.073
#> 7           prof integer  1.00  20.000  44.000  70.500  94.00 4.54e+01  27.509
#> 8   eval_rounded numeric  2.00   4.000   4.000   4.000   5.00 4.00e+00   0.603
#> 9 eval_rounded_r numeric  2.00   3.000   3.000   3.000   5.00 3.00e+00   0.603
#>     n missing
#> 1 463       0
#> 2 463       0
#> 3 463       0
#> 4 463       0
#> 5 463       0
#> 6 463       0
#> 7 463       0
#> 8 463       0
#> 9 463       0

tally(tenure~gender, data = TeachingRatings)
#>       gender
#> tenure female male
#>    no      50   52
#>    yes    145  216

tally(~gender, data = TeachingRatings)
#> gender
#> female   male 
#>    195    268

diffprop(tenure~gender, data = TeachingRatings)
#> diffprop 
#>  -0.0624

TeachingRatings$eval[1] <- NA
TeachingRatings$gender[1] <- NA
head(TeachingRatings$eval)
#> [1]  NA 4.5 3.7 4.3 4.4 4.2

library(sjmisc)

descr(TeachingRatings, eval, beauty)
#> 
#> ## Basic descriptive statistics
#> 
#>     var    type  label   n NA.prc mean   sd   se    md trimmed
#>    eval numeric   eval 462   0.22    4 0.56 0.03  4.00    4.03
#>  beauty numeric beauty 463   0.00    0 0.79 0.04 -0.07   -0.05
#>              range iqr  skew
#>        2.9 (2.1-5) 0.8 -0.46
#>  3.42 (-1.45-1.97) 1.2  0.52
frq(TeachingRatings, gender)
#> 
#> gender <categorical>
#> # total N=463  valid N=462  mean=1.58  sd=0.49
#> 
#> Value  |   N | Raw % | Valid % | Cum. %
#> ---------------------------------------
#> female | 194 | 41.90 |   41.99 |  41.99
#> male   | 268 | 57.88 |   58.01 | 100.00
#> <NA>   |   1 |  0.22 |    <NA> |   <NA>

library(corrr)

TeachingRatings %>% 
  select(eval, beauty, age) %>% 
  correlate()
#> # A tibble: 3 x 4
#>   rowname    eval beauty     age
#>   <chr>     <dbl>  <dbl>   <dbl>
#> 1 eval    NA       0.189 -0.0503
#> 2 beauty   0.189  NA     -0.298 
#> 3 age     -0.0503 -0.298 NA

6.8 Zeilenmittelwerte bilden

Bei Umfragen kommt es häufig vor, dass man Zeilenmittelwerte bildet. Wieso? Man möchte z.B. in einer Mitarbeiterbefragung den “Engagementwert” jedes Beschäftigten wissen (klingt einfach gut). Dazu addiert man die Werte jedes passenden Items auf. Diese Summe teilen Sie durch die Anzahl der Spalten

💡 Zeilenmittelwerte kann man übrigens auch mit Excel bilden …

In R gibt es im Paket sjmisc einen komfortablen Befehl, row_means(). Das Ergebnis sieht man in Tabelle 6.1.

TeachingRatings <- row_means(TeachingRatings, 
                             age, eval, beauty,  # aufzusummierende Spalten 
                             n = .9,
                             var  = "Zeilenmittel")

Der Parameter n gibt den Mindestanteil nicht-fehlender Daten an (hier: 90%), damit ein Zeilenmittelwert berechnet wird. Mit var geben wir einen Namen für neue Spalte des Zeilenmittelwerts an.

Der Befehl row_sums() leistet Ähnliches für Summen- anstatt für Mittelwerte.

6.9 Daten bivariat zusammenfassen

cor(eval ~ beauty, data = TeachingRatings)

TR2 <- dplyr::select(TeachingRatings, eval, beauty, age)
cor(TR2)
#>        eval beauty    age
#> eval      1     NA     NA
#> beauty   NA  1.000 -0.298
#> age      NA -0.298  1.000

ggpairs(TR2)

diffprop(minority ~ gender, data = TeachingRatings)
#> diffprop 
#>   0.0759

6.10 Zeilen filtern

Ist man daran interessiert, nur einen Teil der Fälle (=Zeilen) auszuwerten, so hilft der Befehl filter weiter; filter wird über das Paket tidyverse geladen.

filter(TeachingRatings, gender == "male") -> dozi_maenner

💻 AUFGABE:

• Erstellen Sie eine Tabelle mit festangestellten Dozenten !
• Erstellen Sie eine Tabelle nur mit gut aussehenden Dozenten (der genaue Wert bleibt Ihnen überlassen).

6.11 Spalten auswählen

Manchmal hat meine “breite” Tabelle, also viele Spalten. Da hilft Abspecken, um die Sachlage übersichtlicher zu machen. Sprich: Nur ein paar wichtige Spalten auswählen, die anderen unter den Tisch fallen lassen.

Das kann man mit dem Befehl select (engl. auswählen) erreichen, der über das Paket tidyverse geladen wird:

select(TeachingRatings, eval, beauty) -> TR2

Der Befehl kann noch ein paar Tricks, die man z.B. hier nachlesen kann.

💻 AUFGABE:

• Erstellen Sie eine Tabelle nur mit gender und tenure. Dann wenden Sie tally darauf an.
• Wenden Sie dann die Befehle rowSums und rowMeans auf eine andere von Ihnen erstellten “Mini-Tabelle” an. Speichern Sie das Ergebnis von rowSums als neue Spalte von TeachingRatings.

6.12 Spalten einer Tabelle sortieren

Bestimmt haben Sie schon mal in Excel eine Spalte sortiert, z.B. so, dass die großen Eurowerte ganz oben standen. In R kann man das mit dem Befehl arrange (via tidyverse) erreichen:

arrange(TeachingRatings, -eval) %>% head

Der Befehl %>% head bedeutet nur “UND DANN (das ist das %>%) zeige den Kopf (den Beginn) von dem, was Du gerade gemacht (sortiert) hast”.

💻 AUFGABE:

• Sortieren Sie TeachingRatings nach Schönheit!
• Sortieren Sie TeachingRatings nach Bewertungsergebnis! Sortieren Sie TeachingRatings gleichzeitig nach Schönheit und Bewertungsergebnis! (Tipp: arrange(tabelle, spalte1, spalte2)).

7.1 Syntax von gf_XXX

Die normale Denkweise von mosaic wird verwendet:

gf_diagrammtyp(Y_Achse ~ X_Achse, sonstiges, data = meine_daten).

gf_ steht für ggplot und formula.

Darüber hinaus verkraftet der Befehl noch viele andere Schnörkel, die wir uns hier sparen. Interessierte können googeln… Es ist ein sehr mächtiger Befehl, der sehr ansprechende Diagramme erzeugen kann.

Probieren wir’s!

data(mtcars)
gf_point(eval ~ gender, data = TeachingRatings)

Easy, oder?⁵

Ein anderes Geom:

gf_boxplot(eval ~ gender, data = TeachingRatings)

⚠️ Beachten Sie, dass nur dann mehrere Boxplots gezeichnet werden, wenn auf der X-Achse eine nominal skalierte Variable steht.

Oder mal nur eine Variable (ihre Verteilung) malen:

gf_histogram(~eval, data = TeachingRatings)

💡 Geben wir keine Y-Variable an, nimmt R eigenständig die Häufigkeit pro X-Wert!

7.2 Jittern

Probieren Sie mal diesen Befehl:

gf_point(eval ~ gender, data = TeachingRatings)

Was nicht so schön bei diesem Diagramm ist, ist, dass viele Punkte sich gegenseitig überdecken. Dieses Überdecken bezeichnet man auch als “Overplotting” (hört sich cooler an). Besser wäre es, wenn sich die Punkte nicht überdecken würden, dann würde man besser erkennen, wie viele Punkte wo liegen. Eine einfache Lösung bestünde darin, das Bild etwas zu “schütteln” oder zu “wackeln”, so dass die Punkte etwas verwackelt würden und damit nebeneinander zu liegen kämen. Das kann mit man mit dem Geom jitter (eng. für wackeln) erreichen:

gf_jitter(eval ~ gender, data = TeachingRatings)

Möchte man die Punkte etwas enger haben, so kann man den Parameter width hinzufügen:

gf_jitter(eval ~ gender, data = TeachingRatings, width = .1)

💡 Die Reihenfolge der Parameter in einem R-Befehl ist egal, solange man die Parameter benennt (width, data,…).

💻 AUFGABE:

• Tauschen Sie mal “histogram” mit “density”!
• Erstellen Sie ein Histogramm für beauty!
• Erstellen Sie Boxplots für beauty, vergleichen Sie dabei Männer und Frauen (Tipp: gender steht auf der X-Achse).
• Erstellen Sie Boxplots für eval, vergleichen Sie dabei die überdurchschnittlich schöne mit unterdurchschnittlichen schönen.

7.3 Plot, um Mittelwerte darzustellen

Möchte man nur zwei Mittelwerte darstellen, ist ein Diagramm überflüssig, streng genommen. Schöner ist es, mehr Informationen darzustellen, also z.B. die Rohdaten. Schauen wir uns ein Beispiel aus dem Datensatz tips an:

gf_point(eval ~ gender,
         data = TeachingRatings,
         stat = "summary",
         color = "red", 
         size = 5)

Überprüfen wir mal, ob die Punkte beim Mittelwert liegen:

mean(eval ~ gender, data = TeachingRatings)
#> female   male 
#>   3.90   4.07

😄.

Wir können auch mehrere Gruppen in “Teil-Bildchen” vergleichen, dazu nehmen wir den Operator |; das kann man sich gut merken, wenn man sich vorstellt, dieser vertikale Strich grenzt das linke vom rechten Bild ab:

gf_point(eval ~ gender | minority,
         data = TeachingRatings,
         stat = "summary",
         color = "red", size = 5)

7.4 Wann welches Diagramm?

Ein kurze Übersicht, wann sich welches Diagramm anbietet:

• Mittelwerte vergleichen – Mittelwerte (plus Rohdaten oder Streuungswerte) pro Gruppe darstellen
• Mediane vergleiche – Boxplot
• Verteilung verschiedener Gruppen darstellen – Boxplot (evtl. plus Mittelwert)
• Verteilung einer Gruppe – Dichtediagramm und/oder Histogramm bzw. Balkendiagramm
• Zusammenhang zweier metrischer Variablen – Streudiagramm
• Zusammenhänge nominaler Variablen (Häufigkeiten) – Fliesendiagramm

gf_point(beauty ~ eval, data = TeachingRatings)

gf_point(beauty ~ eval, data = TeachingRatings,
         color = ~ gender,
         shape = ~ gender)

gf_density2d(beauty ~ eval, data = TeachingRatings)  # 2D-Dichte
gf_hex(beauty ~ eval, data = TeachingRatings)  # Kacheldiagramm

tally(minority ~ gender, data = TeachingRatings)
#>         gender
#> minority female male <NA>
#>      no     159  240    0
#>      yes     35   28    1
tally(gender ~ minority, data = TeachingRatings, format = "percent")
#>         minority
#> gender      no   yes
#>   female 39.85 54.69
#>   male   60.15 43.75
#>   <NA>    0.00  1.56

gf_bar(~ minority, data = TeachingRatings) 
gf_bar(~ gender, data = TeachingRatings) 

gf_bar(~ gender, data = TeachingRatings, fill = ~minority, position = "fill") 

7.5 Die Pfeife schlägt zu

Was bedeutet das komische Symbol %>%, welches die beiden Befehle offenbar verkettet? Man nennt es “DIE PFEIFE” (Großbuchstaben machen es erst richtig bedeutsam). Und auf Deutsch heißt dieser Befehl in etwa “UND DANN MACHE…”. Hier verkettet die Pfeife die Beiden Diagrammbefehle, so dass beide Diagramme übereinander gezeichnet werden - ähnlich wie eine Klarsichtfolie, die über ein Bild gelegt wird.

Der Parameter stat = summary führt dazu, dass als Punkte nicht die Rohdaten, sondern eben eine Zusammenfassung (engl. summary) dargestellt wird. In der Voreinstellung ist das der Mittelwert.

Kombinieren wir mal die Rohdaten mit den zugehörigen Mittelwerten in einen Plot:

gf_jitter(eval ~ gender,
          width = .2,
          alpha = .5,  # halb durchsichtig
          data = TeachingRatings) %>% 
gf_point(eval ~ gender,
         data = TeachingRatings,
         stat = "summary",
         color = "red", size = 5)

Diese Sequenz kann grob ins Deutsche übersetzt werden mit:

Hey R,
Erstelle mal die Jitter-Grafik UND DANN  
noch eine Punkt-Grafik für die beiden Mittelwerte (in Rot).

Diese Sequenz ist auch häufig praktisch:

gf_point(eval ~ beauty, data = TeachingRatings,
         color = ~ gender) %>% 
  gf_lm()

7.6 Vertiefung: Ornamente

Hat man einiges an Fingerfertigkeit mit den Diagrammen erlangt, so erwachen Gelüste wie Achsen- oder Titelformatierung oder Farbwahl. Hier einige Hinweise dazu.

gf_point(eval ~ beauty, data = TeachingRatings,
         color = ~ gender) %>% 
  gf_lm() %>% 
  # jetzt die Labels:
  gf_labs(x = "Schönheit",
          y = "Beurteilung des Dozenten",
          title = "Der Zusammenhang von Schönheit und Beurteilung") %>% 
  # jetzt die Farbskala:
  gf_refine(scale_color_viridis_d()) %>% 
  # jetzt die allgemeine Vorlage:
  gf_theme(theme_light())

7.7 Weitere Geome

Hier finden Sie einen Überblick zu Geomen von ggplot, z.B.:

• Boxplot gf_boxplot
• Punkte gf_point
• Linien gf_line
• Histogramm gf_histogram
• …

🔖 Lesen Sie hier weiter, um Ihr Wissen zu vertiefen zu diesem Thema: Vertiefung zur Datenvisualisierung

8.1 Der p-Wert

Ach ja, der p-Wert. Generationen von DozentenStudierenden haben sich wegen ihm oder ob ihm die Haare gerauft. Was war noch mal die Definition des p-Werts? Oder einfacher vielleicht, was will uns der p-Wert sagen?

Der p-Wert gibt an, wie plausibel die Daten unter der getesteten Hypothese sind (der \(H_0\)).

Etwas präziser ausgedrückt:

Der p-Wert gibt die Häufigkeit an, ein Ergebnis, das mindestens so extrem ist, zu bekommen, wenn man den Versuch unendlich oft unter gleichen Bedingungen wiederholen würde.

Gut am p-Wert ist, dass er ein Entscheidungsmaß bietet. Die Gefahr am p-Wert ist, dass man ihn missversteht: Der p-Wert gibt nicht (Sie haben richtig gelesen: nicht) die Wahrscheinlichkeit an, mit der die H0 gilt. Er gibt auch nicht an, wie wahrscheinlich die H1 ist. Er gibt auch nicht an, ob das Ergebnis praktisch bedeutsam ist.

8.2 Simulationsbasierte Inferenz

Simulationsbasierte Inferenz beruht auf einem komputationalen Ansatz: Man führt ein Studie unter kontrollierten Bedingungen mit Hilfe des Computers durch. “Kontrolliert” heißt dabei, dass man eine mit der \(H_0\) kompatible (Daten-)Situation mittels Computer schafft, dann per Zufall auf dieser Grundlage viele Stichproben erstellt und jeweils eine Teststatistik berechnet. Schließlich vergleicht man, wie häufig der beobachte empirische Wert der Teststatistik in diesen nachgestellten (simulierten) Daten ist. Findet sich der echte empirische Wert häufig in den laut \(H_0\) simulierten Daten, so haben wir Grund zur Annahme, dass die Daten mit der \(H_0\) kompatibel sind – auf dieser Basis kann man sich entschließen, die \(H_0\) nicht zu verwerfen. Ist der echte empirische Wert hingegen selten in den simulierten Daten, hat man Grund zu glauben, dass der echte empirische Wert nicht aus der Population laut \(H_0\) stammt. Man verwirft dann die \(H_0\).

Die simulationsbasierte Inferenz hat einige Vorteile:

• es gibt nur allgemeines Prinzip, nicht viele einzelne Verfahren
• es ist voraussetzungsarm; die Normalverteilungsannahme spielt eine geringere Rolle als bei den “klassischen” Verfahren
• das Vorgehen ist einfach und (vergleichsweise) einfach zu verstehen; das geht soweit, dass es Situationen gibt, in denen keine Formeln bekannt sind (auf denen der “klassische Ansatz” der Inferenz aufbaut), man aber mit der Simulation relativ einfach zu einem Ergebnis kommen kann.

Die Nachteile der simulationsbasierten Inferenz sind:

• Man braucht einen Computer, da die Berechnungen von Hand zu umständlich wären
• Bei sehr kleinen Stichproben (\(n < 35\)) sollten die “klassischen” (formelbasierten) Verfahren bevorzugt werden.
• Da der simulationsbasierte Ansatz vergleichsweise neu ist, ist er nicht so bekannt, wie der der klassische Ansatz. Es empfiehlt sich daher, die “Übersetzung” in die klassische (alte) Sprache der Inferenzstatistik zu kennen, das entsprechende Testverfahren (wie t-Test etc.) zu kennen.

8.2.1 Permutationstest

Im Permutationstest wird die Verteilung der Teststatistik unter \(H_0\) simuliert (s. Abbildung 8.1), um den Zusammenhang zweier Variablen zu prüfen. D. h. laut \(H_0\) gibt es gibt keinen Zusammenhang zwischen den zwei Variablen. Untersucht man z.B. ob es einen Zusammenhang im Aufbau einer Webseite (mit Bilder vs. ohne Bildern) und der Verweildauer auf der Webseite gibt, so lautet die \(H_0\), dass sich die Verteilung der Verweildauer in den beiden Gruppen gleich ist.

Hier z.B. für \(H_0: \pi_A-\pi_B=0\):

• Wiederhole oft (z.B. \(10000 \times\)):
  • Mische die \(n_A+n_B\) Beobachtungen.
  • Ordne zufällig \(n_A\) Beobachtungen der ersten Stichprobe zu, die restlichen der zweiten.
  • Berechne die Differenz der zwei Mittelwerte Analoges gilt für andere Teststatistiken, z.B. Anteilsdifferenzen.
• Zeichne ein Histogramm für die Differenz der 10000 Differenzwerte.
• Zähle, wie oft der beobachtete Wert der Teststatistik (oder noch extremere) in der simulierten Verteilung vorkommt.
• Der p-Wert ist der Anteil der zufälligen Teststatistiken (d.h. laut \(H_0\)), die mindestens so groß sind wie der beobachtete Wert.⁶

Abbildung 8.1: Sinnbild zum Permutationstest

8.2.2 Bootstrapping

Beim Bootstrapping wird die Stichprobenverteilung simuliert, z.B. um den Standardfehler \(se\) oder das Konfidenzintervall zu bestimmen (s. Abbildung 8.2).

Hier z. B. für den Mittelwert:

• Wiederhole z.B. \(10000 \times\)
  • Ziehe mit Zurücklegen eine Stichprobe vom Umfang \(n\) aus der Originalstichprobe.
  • Berechne Statistik, z. B. Mittelwert \(\bar{x}\) der Bootstrap-Stichprobe. Analog für andere Statistiken, z.B. Anteil.
• Zeichne das Histogramm der Bootstrap-Verteilung der Statistik.
• Das \(95\,\%\)-Bootstrap-Perzentil-Intervall sind die mittleren \(95\,\%\) der Bootstrap-Verteilung.

Abbildung 8.2: Sinnbild zum Bootstrap

8.3 Auswahl klassischer Testverfahren

Eine kleine Auswahl an Testverfahren wird nachfolgend gezeigt. Diese basieren i. d. R. auf Verteilungsannahmen oder sind asymptotisch und/oder approximativ. Der zuvor gezeigte Permutationstest ist ein exakter Test.

• Binomial-Test: Testet den Anteil \(\pi\) einer (kategorialen) Verteilung.
• Chi-Quadrat-Unabhängigkeitstest: Testet die Unabhängigkeit von zwei kategorialen Variablen.
• t-Test für eine Stichprobe: Testet den Mittelwert \(\mu\) einer Normalverteilung.
• Gepaarter t-Test: Testet die Differenz zweier Merkmale einer Stichprobe unter der Annahme einer Normalverteilung.
• Zweistichproben-t-Test: Testet die Differenz der Mittelwerte zweier Stichproben unter der Annahme einer identischen Normalverteilung.
• Varianzanalyse, ANOVA: Testet die Gleichheit der Mittelwerte zweier oder mehr Stichproben unter der Annahme einer identischen Normalverteilung.
• Korrelationstest: Testet die Korrelation zweier numerischer Variablen einer Stichprobe unter der Annahme einer identischen Normalverteilung.
• Shapiro-Wilk-Test: Testet die Annahme einer Normalverteilung.
• Bartlett’s Test: Testet die Annahme von gleichen Varianzen unter der Annahme einer Normalverteilung.
• Wilcoxon-Test: Testet die Gleicheit von zwei Verteilungen unter der Annahme einer stetigen Verteilung gegen eine Lageverschiebung.
• Kruskal-Wallis Test: Testet die Gleicheit von zwei oder mehr Verteilungen unter der Annahme einer stetigen Verteilung gegen eine Lageverschiebung.

Für eine Mindmap der verschiedenen Verfahren siehe z. B. hier.

Die R-Befehle für die Testverfahren (ob simulationsbasiert oder parametrisch) folgen der normalen Syntax von mosaic:

verfahren(output ~ input, data = meine_daten)

Einige Beispiele:

t.test(eval ~ gender, data = TeachingRatings)
aov(eval ~ gender + native, data = TeachingRatings)
xchisq.test(gender ~ tenure, data = TeachingRatings)

Manchmal ist es sinnvoll, sich das Ergebnis des Verfahrens als Objekt speichern zu lassen und sich dann mit summary() die zentralen Ergebnisse ausgeben zu lassen:

meine_aov <- aov(eval ~ gender + native, data = TeachingRatings)
summary(meine_aov)
#>              Df Sum Sq Mean Sq F value Pr(>F)    
#> gender        1    3.3    3.25   10.98 0.0010 ***
#> native        1    2.8    2.81    9.47 0.0022 ** 
#> Residuals   459  136.1    0.30                   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 1 observation deleted due to missingness

Hinweis: Um den \(\chi^2\)-Wert aus dem Befehl xchisq.test() zu extrahieren, kann man diese Syntax verwenden:

library(mosaic)
chisq_wert <- xchisq.test(smoker ~ day, data = tips)$statistic
chisq_wert

#> X-squared 
#>      25.8

8.4 Die Regression als Schweizer Taschenmesser

💡 Das Schweizer Taschenmesser 🔪 und den Modellierungsverfahren ist die Regressionsanalyse. Man kann sie für viele Zwecke einsetzen.

Weil die Regression so praktisch ist, hier ein Beispiel.

lm(eval ~ beauty, data = TeachingRatings) 
#> 
#> Call:
#> lm(formula = eval ~ beauty, data = TeachingRatings)
#> 
#> Coefficients:
#> (Intercept)       beauty  
#>       3.998        0.133

lm heißt “lineares Modell” - weil man bei der (normalen) Regression eine Gerade in die Punktewolke der Daten legt, um den Trend zu abzuschätzen. Als nächstes gibt man die “Ziel-Variable” (Output) an, hier eval. Dann kommt ein Kringel ~ gefolgt von einer (mehr) Input-Variablen (Prädiktoren, UVs, hier beauty). Schließlich muss noch die Datentabelle erwähnt werden.

Das Ergebnis sagt uns, dass pro Stufe von Beauty die Variable eval um etwa .1 Punkte steigt. Also: Je schöner, desto “besser” sind die Dozenten auch. Immer im Schnitt, versteht sich. (Und wenn die Voraussetzungen erfüllt sind, aber darum kümmern wir uns jetzt nicht.)

Allgemein:

lm(output ~ input, data = meine_daten)

Easy, oder?

Man kann auch mehrere Prädiktoren anführen:

lm(eval ~ beauty + gender, data = TeachingRatings)
#> 
#> Call:
#> lm(formula = eval ~ beauty + gender, data = TeachingRatings)
#> 
#> Coefficients:
#> (Intercept)       beauty   gendermale  
#>       3.882        0.148        0.200

Möchte man ein ausführliches Ergebnis bekommen, so verlangt man von R eine Zusammenfassung (summary) des lm-Befehl, und zwar mit dem Befehl summary. Den Befehl summary kann man mit dem Und-danach-Befehl (%>%) an den lm-Befehl anschließen:

lm(eval ~ beauty + gender, data = TeachingRatings) %>% summary()
#> 
#> Call:
#> lm(formula = eval ~ beauty + gender, data = TeachingRatings)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -1.8721 -0.3700  0.0341  0.4010  1.0322 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   3.8819     0.0388  100.11  < 2e-16 ***
#> beauty        0.1484     0.0320    4.64  4.5e-06 ***
#> gendermale    0.1997     0.0511    3.91  0.00011 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.538 on 459 degrees of freedom
#>   (1 observation deleted due to missingness)
#> Multiple R-squared:  0.0667, Adjusted R-squared:  0.0626 
#> F-statistic: 16.4 on 2 and 459 DF,  p-value: 1.32e-07

Dazu werden die durch + getrennt. Pro Prädiktor wird die Steigung der Regressionsgeraden angegeben. Man kann auch nominale Prädiktoren reinfüttern. Das macht die Regression so praktisch.

In diesem Fall sehen wir, dass Schönheit einen positiven Koeffizienten aufweist, d.h. die Regressionsgerade steigt: Für jeden Punkt Schönheit steigt die (mittlere) Bewertung um etwa 0.15 Punkte. Für Geschlecht gilt, dass genderFemale (die Frauen) um etwa 0.20 schlechter (wegen dem Minuszeichen) in der Beurteilung eingeschätzt werden.

gf_point(eval ~ beauty, 
         data = TeachingRatings,
         color = ~gender) %>% 
  gf_lm()

lm(eval ~ beauty + gender + gender:beauty, data = TeachingRatings) %>% 
  summary()
#> 
#> Call:
#> lm(formula = eval ~ beauty + gender + gender:beauty, data = TeachingRatings)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -1.8382 -0.3746  0.0415  0.3999  1.0676 
#> 
#> Coefficients:
#>                   Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)         3.8889     0.0389  100.00  < 2e-16 ***
#> beauty              0.0871     0.0471    1.85  0.06504 .  
#> gendermale          0.1970     0.0510    3.86  0.00013 ***
#> beauty:gendermale   0.1132     0.0640    1.77  0.07773 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.536 on 458 degrees of freedom
#>   (1 observation deleted due to missingness)
#> Multiple R-squared:  0.073,  Adjusted R-squared:  0.0669 
#> F-statistic:   12 on 3 and 458 DF,  p-value: 1.37e-07

mein_lm <- lm(eval ~ beauty, data = TeachingRatings)

predict(mein_lm, data.frame(beauty = 1.5))
#>   1 
#> 4.2

vohersagefehler <- residuals(mein_lm)

RMSE <- sqrt(mean(vohersagefehler^2))
RMSE  
#> [1] 0.545

8.5 Effektstärkemaße

Der p-Wert ist ein Maß für die Plausibilität einer Hypothese; er ist kein Maß für die Stärke eines Effekts. Grund dafür ist u.a., dass die Stichprobengröße in die Berechnung des p-Werts einfließt. Große Stichproben werden daher (unter sonst gleichen Umständen) schneller signifikant als kleine Stichproben (und sehr große viel schneller…).

Um die Stärke von Effekten zu bemessen, gibt man den Wert der Teststatistik an (oder eine aus ihr abgeleitete Größe); diese bezeichnet man als Effektstärkemaße. Vergleicht man z.B. zwei Mittelwerte ist diese Mittelwertsdifferenz die Effektstärke; teilt man diesen Wert noch durch die mittlere Standardabweichung, so erhält man eine häufig verwendete Kenngröße, Cohens d. Berechnen wir Cohens d für den Mittelwertsunterschied der Beurteilung abhängig vom Geschlecht:

library(lsr)
cohensD(eval ~ gender, data = TeachingRatings)
#> [1] 0.309

Bei einer Varianzanalyse ist das Verhältnis von erklärter zu nicht erklärter Varianz oder zu Gesamtvarianz die typische Effektgröße. Diese Kenngröße wird als \(\eta^2\) (Eta-Quadrat) bezeichnet oder auch als \(R^2\). Ebenfalls im Paket lsr gibt es einen Befehl für diese Effektgröße:

meine_aov <- aov(eval ~ gender, data = TeachingRatings)

etaSquared(meine_aov)
#>        eta.sq eta.sq.part
#> gender 0.0229      0.0229

Effektstärken für Anteilsdifferenzen werden am besten als Anteilsdifferenzen wiedergegeben, z.B. im Rahmen von \(\chi^2\)-Tests. Bei \(\chi^2\)-Tests, bei denen mindestens eine der Variablen mehr als zwei Stufen aufweist, kann auf Cramers V ausgewichen werden:

cramersV(TeachingRatings$native, TeachingRatings$tenure)
#> [1] 0.124

Bei (linearen) Regressionsanalysen ist \(R^2\) ein Effektstärkemaß; es quantifiziert den Anteil der durch die Prädiktoren erklärte Varianz. Die Wichtigkeit eines Prädiktors kann man z.B. am unstandardisierten Einflussgewicht (Betakoeffizient) ablesen. Gängig ist auch, den Zuwachs an \(R^2\) zu berichten, wenn man einen Prädiktor zu einem Regressionsmodell hinzufügt.

Bei einer Korrelationsanalyse ist die Sache einfach: Der (absolute) Wert von \(r\) ist das Effektstärkemaß.

9.1 Tabellen und Diagramme

Daten kommunizieren heißt praktisch zumeist, Tabellen oder Diagramme zu erstellen. Meist gibt es dazu Richtlinien von Seiten irgendeiner (selbsternannten) Autorität wie Dozenten oder Fachgesellschaften. Zum Beispiel hat die APA ein umfangreiches Manual zum Thema Manuskriptgestaltung publiziert; die deutsche Fachgesellschaft der Psychologie entsprechend. Recherchieren Sie mal, wie in Ihren Richtlinien Tabellen und Diagramme zu erstellen sind (oder fragen Sie Ihren Gutachter). Es gibt einige R-Pakete, die helfen, Tabellen zu erstellen, etwa apaTables, sjPlot, xtable oder kableExtra.

9.2 Für Fortgeschrittene: RMarkdown

Wäre das nicht cool: Jegliches Formatieren wird automatisch übernommen und sogar so, das es schick aussieht? Außerdem wird Ihr R-Code und dessen Ergebnisse (Tabellen und Diagramme oder reine Zahlen) automatisch in Ihr Dokument übernommen. Keine Copy-Paste-Fehler mehr. Keine händisches Aktualisieren, weil Sie Daten oder die vorhergehende Analyse geändert haben. Hört sich gut an? Ist es auch. Probieren Sie mal RMarkdown aus.